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A numerical simulation of kinetic plasma turbulence is performed to assess the applicability of critical balance 
to kinetic, dissipation scale turbulence. The analysis is performed in the frequency domain to obviate compli- 
cations inherent in performing a local analysis of turbulence. A theoretical model of dissipation scale critical 
balance is constructed and compared to simulation results, and excellent agreement is found. This result 
constitutes the first evidence of critical balance in a kinetic turbulence simulation and provides evidence of an 
anisotropic turbulence cascade extending into the dissipation range. We also perform an Eulerian frequency 
analysis of the simulation data and compare it to the results of a previous study of magnetohydrodynamic 
turbulence simulations. 



I. INTRODUCTION 

Plasma turbulence is ubiquitous in a wide range of 
space and astrophysical environments, playing a funda- 
mental role in transferring energy from the large scales at 
which the turbulence is driven to the small scales at which 
the turbulence is dissipated. Developing a detailed un- 
derstanding of plasma turbulence is one of the key goals 
of the space physics and astrophysics communities. 

One of the central tenets of plasma turbulence is the 
concept of critical balance. Critical balance is the sup- 
position that the time scale associated with linear fluc- 
tuations of Alfven waves, to = k^VA, is of the same order 
of magnitude as time scale associated with the nonlinear 
cascade of energy, ui n i ~ k±v±, where va is the Alfven 
speed, v± is the perpendicular fluctuation velocity, and 
parallel and perpendicular are defined with respect to the 
direction of the local mean magnetic field^— . Critical 
balance leads to the prediction of an anisotropic cascade 
of energy in wavevector space, where magnetic energy 
cascades at different rates parallel and perpendicular to 
the local mean magnetic field. 

Although the original predictions of critical balance 
pertain only to a cascade of Alfven waves in the magne- 
tohydrodynamic (MHD) limit of the inertial range, the 
theory can be extended to scales smaller than the ion gy- 
roradius, at which wave-particle damping and collisions 
become important^—. The latter region is often referred 
to as the dissipation range of plasma turbulence, where it 
is proposed that the cascade of Alfven waves transitions 
to a cascade of kinetic Alfven waves (KAW)£~— . 

The anisotropic scaling of the magnetic field energy 
has been observed in the inertial range portion of the 
solar wind through the use of wavelets or second-order 
structure functions to discern the local mean magnetic 
field direction, e.g.j2r— . Horbury et al.— demonstrated 
that critical balance fits the solar wind observations well 
in the directions parallel and perpendicular to the local 
mean magnetic field, and Forman et al.— went further, 
demonstrating that solar wind observations follow the 



a ' Electronic mail: jason-tenbarge@uiowa.edu 



predictions of critical balance for all angles between par- 
allel and perpendicular. 

Critical balance and its physical consequences have 
also been tested and verified in many numerical turbu- 
lence simulations. It has been tested extensively in MHD 
simulations, e.g.^^—, and at smaller, dissipation range 
scales in electron MHD (EMHD) simulations, e.g.,— 
However, all of the previous studies have been performed 
with fluid codes that cannot capture wave-particle inter- 
actions nor accurately model collisional effects, both of 
which play important roles in dissipating turbulence at 
scales below the ion gyroradius. 

To evaluate whether critical balance persists in the dis- 
sipation range, we consider here a detailed study of a 
numerical simulation of dissipation scale turbulence per- 
formed using AstroGK, the Astrophysical Gyrokinetics 
Code, developed to study kinetic turbulence in astro- 
physical environments. Rather than examining the en- 
ergy distribution in wavenumber space, we perform our 
analysis in the frequency domain to obviate some of the 
difficulties inherent in performing a local analysis of tur- 
bulence, which are discussed in ^VII The frequency is 
used as a proxy for the parallel wavenumber to deter- 
mine whether or not an anisotropic cascade, consistent 
with that predicted for critically balanced kinetic Alfven 
wave turbulence, exists in the dissipation range. An Eu- 
lerian frequency analysis of the AstroGK simulation will 
also be compared to a similar study by Dmitruk and 
Matthaeus— performed using a MHD simulation. 

II. SIMULATION PARAMETERS 

A detailed description of AstroGK and the results of 
linear and nonlinear benchmarks are presented in Nu- 
mata et al.—, so we provide here only a brief overview. 

AstroGK is an Eulerian slab code with triply periodic 
boundary conditions that solves the electromagnetic gy- 
roaveraged Vlasov-Maxwell five dimensional system of 
equations. It solves the gyrokinetic equation and gy- 
roaveraged Maxwell's equations for the perturbed gy- 
roaveraged distribution function, h s (x, y, z, A, e), for each 
species s (protons and electrons), the parallel vector po- 
tential, A z , and perturbed parallel magnetic field, SB Z , 
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and the scalar potential, ^ i 22 ' 23 . The simulation domain 
is elongated along the direction of the equilibrium mag- 
netic field Bo = BqZ. Velocity space coordinates are re- 
lated to the energy, e = v 2 /2, and pitch angle, A = v\/v 2 . 
The equilibrium distribution for both species is treated as 
Maxwellian, and a realistic mass ratio, m p /m e = 1836, is 
used. The x — y plane is treated pscduospectrally, and an 
upwinded finite-diffcrencing approach is employed for the 
z-direction. Velocity space is evaluated following Gaus- 
sian quadrature rules. Linear terms are evolved implic- 
itly in time, while nonlinear terms are evolved explicitly 
by a third-order Adams-Bashforth method. Collisions 
are treated using a fully conservative, linearized, and gy- 
roaveraged collision operato r 24 ' 25 . 

To represent average solar wind conditions observed 
at ~ 1 AU, we choose plasma parameters /3j = 1 and 
Ti = T e , where fa = v 2 Jv\, va = B / y/Awn^rnl, and 
Vu = y // 2Ti/rrii is the ion thermal speed. We output the 
electromagnetic field information at fixed time intervals. 
Since this diagnostic is data intensive, we choose to per- 
form the simulation on a numerical grid smaller than the 
largest simulations performed with AstroGK 2 ^. Specifi- 
cally, we use a numerical grid of (n x , ny,n z ,n e ,n\,n s ) = 
(32, 32, 64, 16, 16, 2), where n e , n\, and n s are the number 
of grid points in energy, pitch angle, and species respec- 
tively. The spatial extent of the domain is (L x , L yi L z ) — 
2n(pi, pi, ao), where pi = Vti/Cli is the ion gyroradius, 
f2j = eBo/rriiC is the ion gyrofrequency, and a deter- 
mines the domain elongation and is chosen by assum- 
ing a critically balanced Alfvenic inertial range: 6 = 
Pi/a® = (kiPi) 1 / 3 {kj_oPi) 2 ^' > , where ki is the wavenumber 
corresponding to the outer scale of the turbulent cascade 
(much larger than our simulation domain) at which the 
turbulent system we are modelling is physically driven, 
and fcj_o is the perpendicular wavenumber corresponding 
to the simulation domain size, the largest resolved length 
scale in the simulation (sub-script is used throughout 
to indicate domain-scale quantities). Based on in situ 
measurements, we choose kiPi = 10~ 4 for the solar wind 
at 1 AU&. Assuming this value for fc, implies 5 ~ 1/20, 
yielding a simulation domain with the z— direction elon- 
gated by a factor of 1/8 — 20. 

Normalization of the time t = t/n and frequency 
u> = LOT t uses the parallel thermal time, r t = ao/vu = 
&o/(f aVWi where = 1. The normalized linear ki- 
netic Alfven wave frequency at the largest scale of our 
simulation is luq = 1.1366, determined from a linear gy- 
rokinctic dispersion relation solver—. The corresponding 
normalized domain-scale time is given by 
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Electromagnetic field data is output every At = 0.02r t , 
resulting in a Nyquist frequency of Cjjq q ~ 160. 

We drive our simulation with an oscillating Langevin 
antenna driving the parallel vector potential. Details re- 
garding the driving antenna are provided in TenBargc 
ct al. (2011)21, so here we specify only the antenna pa- 



rameters. We drive the largest four independent modes 
of our simulation domain, (k x ,k y ,k z /S)pi = (±1,0, ±1) 
and (0,±1,±1), with a frequency uj a ~ ujq and decorre- 
lation rate "f a ~ 0.7wo. The antenna amplitude, A\\ , is 
chosen to satisfy critical balance at the domain-scale, 
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where C2 — 1 is a Kolmogorov constant and an analyt- 
ical estimate of the linear Alfven /kinetic Alfven wave 
frequency normalized to kuVA is given by 
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Plots of the perpendicular magnetic energy spectrum 
at the beginning (red) and end (solid black) of the anal- 
ysis included herein can be found in Figure [U The anal- 
ysis begins after the cascade has become well-developed, 
which corresponds to t = 0.87to, and the analysis ends 
at t = 7.44to. The vertical dotted line at k±pi = 10 cor- 
responds to the maximum fully resolved perpendicular 
mode of the simulation. Beyond that point, hypercol- 
lisionality acts to remove energy from the system and 
steepens the spectrum. Due to the hypercollisionality, 
the domain of validity is limited to k±pi < 8. The satu- 
rated spectra have a spectral index ~ —2.8. This spec- 
tral index agrees with larger scale numerical simulations 
with the same plasma parameters 2 ^. This demonstrates 
that, although the dynamic range of the simulation un- 
der consideration is rather limited, the overall behaviour 
is consistent with a larger scale simulation. 



III. CRITICAL BALANCE THEORY 

The strength of turbulence can be characterized by the 
ratio of the nonlinear cascade rate or perpendicular eddy 
turn-around time, 0J n i ~ fe_i_u_i_, to the linear frequency, 
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When the turbulent fluctuations exist for several 

turn-around times at a given scale before their energy is 
cascaded to smaller scales — many wave- wave "collisions" 
are required to cascade their energ y 28 ! 29 . This situa- 
tion, known as weak turbulence^, generates a cascade of 
energy in only the perpendicular directio n 31 ' 32 . There- 
fore, x grows with increasing perpendicular wavenumber, 
strengthening the nonlinear interactions. Once % ~ 1, 
the timescales of the nonlinear and linear processes be- 
come equal, and the turbulence is said to be critically 
balanced 1 ^. 

It is interesting to consider also the over-strong case, 
^ » 1. In this case, the nonlinear frequency is larger 
than the linear frequency, so two interacting Alfven wave 
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FIG. 1: (Color online) Perpendicular magnetic energy 
spectra at the beginning (red) and end (solid black) of 
the frequency analysis compared to a larger, higher 
resolution simulation with similar plasma parameters 
(blue). The dashed line corresponds to a spectral slope 
of —2.8, and the vertical dotted line corresponds to the 
maximum fully resolved perpendicular scale. 



packets can undergo multiple cascades to smaller scales in 
a single linear crossing time. Therefore, the Alfven wave 
packets are expected to be rapidly cascaded in the par- 
allel direction until they restore the condition of critical 
balance, x — 1- Note that fluctuations with k\\ ~ 0, and 
therefore with u> ~ 0, are naturally generated by three- 
wave interactions of the Alfvenic turbulence 1 -^. This 
regime of over-strong turbulence is dominated by un- 
corrected fluctuations, because fluctuations of scale k± 
are decorrelated over parallel distances greater then k\\x- 
Thus, the energy in this over-strong region of wavenum- 
ber space is expected to be roughly constant, as observed 
in simulations of MHD turbulence^. 

The assumptions of critical balance and constant 
energy cascade rate lead to a predicted wavenumber 
anisotropy scaling of 
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which can be combined into a single equation of the form 
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where £ = 1/3 gives the standard dissipation range scal- 
ing derived from fluid theories assuming critical balance 
holds in the dissipation ranged—. The form of Equa- 
tion (JBJ has been chosen so that the asymptotic limits 
are continuously connected and it is well-behaved in its 
domain. Note that substituting Equation (jSJ) into Equa- 
tion (|3J) yields the linear Alfvcn/KAW frequency u in 
terms of only kj_ . 



A schematic diagram depicting the expected popula- 
tion of energy in critically balanced Alfvenic turbulence 
in terms of ui and k± is presented in Figure The dis- 
persion relation defines the critical balance boundary in 
the k_\_-uj plane and corresponds to the solid line in Figure 
[21 Critical balance suggests the turbulent energy will fill 
the region below the critical balance boundary (shaded), 
as seen in incrtial range MHD simulations of turbulence, 
e _g_ ( 3 i i5- i i8 ) an j m dissipation range electron MHD simu- 
lations, e.g.j2ii&. 
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FIG. 2: Schematic diagram depicting the population of 

energy in the k±-u plane assuming critical balance 
holds. The solid lines correspond to the critical balance 
boundary. The isotropic energy injection scale is 

k iPi = 10" 4 . 

To obtain a more realistic prediction of critical balance, 
we construct a physical model of the expected energy 
distribution. Previous models of critical balance, e.g.j 1 ^, 
assume a form for the energy distribution in fe_i_-fcii space 

of E(k h k x ) ~ fc; 10/3 /(fc||£ 1/3 A-l /3 ), where f(\u\) ~ 1 
for u < 1 and is negligibly small for u^$> 1. 

We assume a similar functional form in k±-u> space, 



E{u,k x ) = 



E, 



(k±Pi)~ 



1 



l + (fc±p 2 ) 2 



1 + (uj/uicb) 



(7) 



where Eq is an overall energy normalization, 7 and e are 
the desired inertial and dissipation range energy scal- 
ings, and w c & = fcjjaow/v 7 /^- The functional form for 
the k± portion of Equation ([7]) is chosen empirically so 
that the frequency integrated energy reduces to a one- 
dimensional perpendicular energy spectrum that scales 
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as EB x (k±) = f E(u>,k±) (ko oc kj 7 in the inertial 
range and £f[ e in the dissipation range. We take 7 = 5/3 
and e = 2.8, consistent with the dissipation range solar 
wind2£ and large scale kinetic numerical simulations^. 
The functional form of the frequency has been chosen to 
have a flat energy spectrum up to ~ uj c b followed by an 
8-th order roll-off at higher frequencies. 

The logarithm of Equation ([7]) is plotted in Figure 
I3al Over plotted in Figure [3a| is the linear frequency, 
which corresponds to the critical balance boundary. To 
make the boundary more clear, normalizing the energy at 
each k± by the zero frequency energy, E{u, k±)/ E(0, k±), 
yields Figure [3b] 

Since we have now constructed a more realistic nu- 
merical model of the critical balance energy distribution, 
given by Equation ([?])■ we can examine the fraction of 
the turbulent energy falling below the critical balance 
boundary (in frequency) relative to the total energy at 
each k±, 

E f . E(u, fci) duj 
Etot J E(uj,k^) did 

This fraction is plotted as the solid line in FigureS] Up to 
the point where finite box size limitations become impor- 
tant, ~ 90% of the energy lies below the critical balance 
boundary. 

IV. CRITICAL BALANCE SIMULATION 

We may now compare the energy distribution in fc^- 
u> space from the theoretical model for critical balance, 
given by Equation (J7J , with that determined through fre- 
quency analysis of our kinetic Alfven wave turbulence 
simulation using AstroGK. 

In Figure [5a] is plotted the logarithm of the perpen- 
dicular magnetic energy. The figure was constructed by 
integrating the 3D AstroGK data in k z and annular sec- 
tions in the k x -k y plane. The over-plotted curves (black) 
correspond to critical balance with £ in Equation 
set to (dashed), 1/3 (solid), and 1 (dot-dashed), where 
£ = 1/3 is the conventional prediction for critically bal- 
anced kinetic Alfven wave turbulence in the dissipation 
range. 

Below the £ = 1/3 (solid black) critical balance bound- 
ary in Figure I5a[ we find excellent qualitative agreement 
with the theoretical prediction represented in Figure [3a] 
Although the agreement is slightly poorer above this 
boundary, the fraction of the energy in this region is 
very small. The energy fraction falling below each of the 
curves in Figure [5a] is plotted in Figure |U The upward 
trend beginning at k±pt ~ 8 is due to hypercollisionality 
and finite box size limitations, and the poor agreement 
at small k± is due to the effect of driving. Away from 
these limiting values, the difference between the theory 
and the simulation for £ = 1/3 is within approximately 
10%. 



In Figure [5b] is plotted the perpendicular magnetic en- 
ergy at each k± normalized to the zero frequency energy, 
EB ± {u,k±)/EB ± (0,k±). The figure was constructed in 
the same manner as Figure I5al and the curves (white) 
correspond to the same values of £. As demonstrated 
with the theoretical model, plotting the energy distribu- 
tion in uj-kj_ space normalized in this fashion highlights 
the critical balance boundary, which closely follows the 
standard dissipation range prediction of £ = 1/3 (solid 
white). Again, the qualitative agreement with the pre- 
diction of critical balance is excellent. The non-smooth 
distribution of energy below the critical balance line in 
Figure I5bl is primarily due to the discrete nature of this 
moderate resolution simulation. 

Another method of visualizing the perpendicular mag- 
netic energy distribution is to perform cuts along Gj at 
fixed fcj_, as presented in Figure [B] Since the energy is 
plotted linearly, the area under each curve corresponds 
to the turbulent energy. Each panel of the figure repre- 
sents a different k± slice, and the vertical dashed lines 
indicate the £ = 1/3 critical balance boundary. Panel 
a is influenced by the driving at k^pi = 1, but it still 
displays similar qualitative behaviour to cuts at higher 
k±pi. The general trend can be seen to be an approxi- 
mately flat energy distribution up to a frequency some- 
what less than the critical balance boundary, followed by 
a steep roll-off. The majority of the energy in all cases 
is contained within the critical balance boundary. This 
plot makes clear an important characteristic of the turbu- 
lence: no energy significantly in excess of that predicted 
by the critical balance model is seen either at uj = or 
at frequencies above critical balance. 



V. EULERIAIM FREQUENCY ANALYSIS 

Eulerian frequency spectra are constructed by placing 
an array of "probes" across a single spatial x-y plane in 
the middle of the simulation domain. We use an array of 
64 x 64 probes across the plane to record a time series of 
the fluctuating magnetic field components. A schematic 
of the probe distribution can be seen in Figure [7J 

The evolution in time of the three magnetic field com- 
ponents at a single probe location over the same time 
interval used for the preceding critical balance study are 
illustrated in the upper panel of Figure[8] For illustrative 
purposes, over-plotted in the figure is cos (t/ro + 0.5). 
Clearly, the B x (blue) component is dominated by <1> ~ 1 
fluctuations, which corresponds to the driving frequency 
and the lowest linear mode of the system. 

In the lower panel of Figure [8] are plotted the Fourier 
response of the three magnetic field components averaged 
over the full probe array. All three magnetic field com- 
ponents are dominated by spectral peaks at uj ~ and 
±1, with a clear band-gap between the values. As noted 
above, Cj ~ ±1 corresponds to the driving frequency and 
the lowest linear mode of the system. As such, this spec- 
tral component should contain the most power due to 
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FIG. 3: (Color online) (a) log [Eb± (w, k±)] and (b) EB ± {u,k±)/EB ± (0,k±), both computed based upon theoretical 
expectations. The white curve in both panels corresponds to the linear dispersion relation of kinetic Alfven waves 

with f = 1/3. 
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FIG. 4: The fraction of energy beneath the critical 
balance boundary at each k±. The solid line is based 
upon theoretical considerations of critical balance, while 
the other three curves are calculated from the AstroGK 
simulation with different values of £. 



the driving, responsible for generating the forward cas- 
cade of energy to higher frequencies. We conjecture that 
the peak at Cj = exhibits significant energy because this 
mode is responsible for nonlinear scattering in three- wave 
interactions of turbulence and is self-consistently gener- 
ated via the nonlinear interactio n 30 ' 31 ! 34 . 



Figure [S] presents the Eulcrian frequency spectrum 
of the perpendicular magnetic energy averaged over all 
probes. The vertical dotted lines indicate the minimum 
and maximum linear frequencies for the minimum and 
maximum resolved length scales in the simulation do- 
main, as given by Equation ©. The frequency range 
around a) ~ 1 is dominated by the antenna driving, so 
we fit from Cj = 2 to 14.4 to obtain a spectral index of 
~ — 3.2 for the Eulerian frequency spectrum. 

Excluding the the driving dominated portion of the 
spectrum around Cj ~ 1 , the energy in the low frequency 
range below Cj^ in is approximately constant, consistent 
with the predictions of critical balance outlined in £|IIII 
It is important to note that very little energy resides in 
these low frequency modes: the total integrated energy 
in the turbulence is the area under a linear plot of the 
frequency spectrum, so this low- frequency range corre- 
sponds to very little integrated area but is exaggerated 
in the logarithmically plotted Figure [9] 



VI. DISCUSSION 

The model for critical balance constructed in i jllll is 
based upon the theoretical expectations for strong tur- 
bulence, the foundation for which was established by Gol- 
dreich and Sridhar— : our model quantifies the qualitative 
arguments given therein and extends the notion of crit- 
ical balance into the dissipation range, permitting arbi- 
trary perpendicular wavevector scaling to agree with the 
results of large-scale kinetic numerical simulations^ and 
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FIG. 5: (Color online) (a) log [Eb ± (w, k±)] and (b) Es ± (uJ,k±)/EB ± (0,k±) from the AstroGK simulation. The 
curves correspond to critical balance with £ = 0, 1/3, and 1 in ascending order. 
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FIG. 6: Perpendicular magnetic energy normalized to 
the zero frequency energy at each k± versus frequency 
at fixed k±. The vertical lines represent the £ = 1/3 
critical balance boundary. Panels a) - e) step through 
kj_ bins. 
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FIG. 7: (Color online) Eb^ at a single z plane in the 
middle of the simulation domain. Overlaid are the 
positions of the "probes" used to record temporal 
fluctuations. 



solar wind observations^ of the dissipation range. 

The theoretical model of critical balance is compared 
to our numerical simulations in i jlVI Considering that the 
concept of critical balance is based upon order of mag- 
nitude scaling arguments, the agreement between theory 
and simulation presented in Figure I5bl is striking. It is 
important to note that the nature of the energy input into 
our simulation domain generates perpendicular magnetic 



energy fluctuations that resemble linear Alfvcn/KAWs. 
The turbulent fluctuations, however, are only driven at 
the smallest resolved wave numbers of our simulation do- 
main, with the frequency range of the driving energy cen- 
tered at Cj = Qjq. Therefore, all of the energy at larger 
wave numbers and frequencies away from ljq results from 
self-consistent, nonlinear turbulent interactions. 

The spectral anisotropy in fluid simulations is typically 
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FIG. 8: (Color online) In the upper panel is presented 
the temporal evolution at a single probe location of the 
three magnetic field components: B x (blue), B y (red), 

and B z (green). The dotted black line is 
cos (t/ro + 0.5). In the lower panel are presented the 
Fourier responses of the three magnetic field 
components averaged over the full probe array. 



determined via two-point, second-order structure func- 
tion analyses, e.g.,— ii i 15 i 35 . This approach allows one to 
identify the wavevector parallel to the local mean mag- 
netic field at a given scale. A standard Fourier analy- 
sis provides only wavevectors parallel and perpendicu- 
lar to the global mean magnetic field, i.e., k z and &r, 
rather than k\\ and k±£^&. This effect is a consequence 
of the inherent averaging of the Fourier transform. Tur- 
bulent fluctuations relevant to solar wind turbulence in 
the dissipation range are characterized by SB±/B <C 1 
and ku/kx < 1, so the Fourier transformed wavevector 
components are given, to lowest order, by kn ~ k± and 
k z ~ 9kj_, where 9 ~ SB±/B^. Thus, a standard Fourier 
approach to analyzing turbulent simulations is sensitive 
only to the spectrum in k±^ and the amplitude of the 
fluctuating field, rather than ku spectrum along the local 
mean magnetic field. 

The structure function approach works well in MHD 
inertial range turbulence simulations and undamped dis- 
sipation range turbulence simulations, wherein the spec- 
tral slope is relatively shallow. The perpendicular spectra 
in kinetic simulations of dissipation range turbulence^ 
and in solar wind observations^ have spectral indices 
around —2.8, and the parallel spectrum is expected to be 
steeper yet^i^S. Two-point, second-order structure func- 
tions cannot recover spectral indices steeper than —3^2, 
and thus structure functions cannot be used to analyze 
the parallel spectrum of our dissipation range simula- 
tions. 

The analysis performed herein obviates the limitations 
of structure function and spatial Fourier analyses by us- 




FIG. 9: The perpendicular magnetic energy spectrum 
as a function of frequency averaged over the full probe 
array. The dashed line corresponds to a spectral index 
of —3.2. The vertical dotted lines correspond to the 
expected minimum and maximum frequencies from 
linear theory. 



ing the frequency w as a proxy for fcii . This approach 
is motivated by the fact that the linear frequency for 
the Alfven and kinetic Alfven wave is linearly propor- 
tional to parallel wavenumber, lj on fen. This relation 
can be seen clearly by rewriting Equation ([3]) in the form 
ui(k±, fcii) = u](k±)k\\VA- The observation of a critical 
balance scaling in the turbulent energy distribution in 
k±_-u space for our kinetic simulation of dissipation range 
turbulence implies that the spectral anisotropy observed 
in MHD inertial range simulation s 3 ' 15 " — and electron 
MHD dissipation range simulations^^ persists even in 
a turbulent kinetic plasma. Since the conventional crit- 
ical balance boundary given by £ = 1/3 is a good fit to 
the turbulent energy distribution from our simulation, we 
conclude that the turbulent spectrum of kinetic Alfven 
waves is well described by an anisotropic, critically bal- 
anced energy distribution in wavevector space given by 
< £ < 1/3, even for a fully nonlinear and collisionlessly 
damped kinetic simulation. 

The results of ^IVI and |V] demonstrate that the the- 
ory of critical balance extends into the dissipation range 
of turbulence. Since critical balance is fundamentally a 
balancing of the linear and nonlinear processes in plasma 
turbulence, our results suggest that linear wave theory 
plays an important role, even in a strongly nonlinear tur- 
bulent plasma. 

Dmitruk and Matthaeus— (DM) performed a similar 
Eulcrian frequency analysis of driven MHD turbulence, 
but their conclusions about the nature of the turbulence 
differ dramatically from our findings. In particular, they 
claim to find an excess accumulation of energy extend- 
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ing from their lowest linear mode to zero frequency and 
a very steep spectrum above their lowest linear mode. 
They conclude that linear waves play little role in MHD 
plasma turbulence. We believe the conclusions of DM 
are significantly biased by the set-up of their numerical 
simulations, for a number of reasons detailed below. 

First, DM drive their simulations isotropically (with 
k±_ ~ feii for the driven modes) with a fixed ampli- 
tude, yielding a turbulent fluctuation amplitude SB ~ 1, 
for a range of values of the mean field strength Bq = 
0, 1, 2, 8, 16. The nonlinearity parameter for Alfvenic tur- 
bulence may be expressed as x = k±5B/k\\Bo, so the only 
case that yields strong MHD turbulence with \ ~ 1 is the 
B = 1 case; all other cases correspond to simulations of 
weak turbulence, with increasingly small nonlinearity pa- 
rameters, x = 1/2, 1/8, 1/16. That the turbulence is in- 
deed weak is supported by the very steep spectral indices 
of their frequency spectra. 

The nonlinear frequency in the weak turbulence 
regime 3 -^ is given by u) n i ~ x 2 k\\VA , which suggests that 
one may expect to sec a signature in the frequency spec- 
trum corresponding to this very low nonlinear frequency, 
as indeed observed by DM. In addition, DM report that 
their simulations require "tens of nonlinear times" to 
reach a saturated steady state, although they do not de- 
fine what they mean by a "nonlinear time." If nonlinear 
time is taken to mean the Alfven crossing time, their re- 
sult is consistent with the expectation that it will take 
approximately one full nonlinear timescale, t„; ~ ta/x , 
or many Alfven crossing times, ta = ^/k\\VA, to reach 
the steady state of the weak turbulent cascade. Note 
that other simulations of driven, strong MHD turbulence 
report the establishment of a steady state within one to 
a few Alfven times, as expected for x ~ 1— i 15 ' 19 ' 39 ' 40 . 
Therefore, we speculate that the excess energy at low 
frequency observed in the weak turbulence simulations 
of DM may be attributed cither to the response of the 
plasma at the low nonlinear frequency or to the inclusion, 
in the interval used for the frequency analysis, of the long 
transient evolution of the turbulence as it approaches a 
steady state. 

Second, in weak turbulence, the resonant match- 
ing conditions for the nonlinear term in wavevec- 
tor and frequenc y 30 ' 32 for the dominant three wave 
interaction o 2 ' 34 ' 41 " — imply a crucial role for modes with 
fc|| = in the nonlinear transfer of energy. Accord- 
ing to the linear dispersion relation for Alfven waves, 
u> = ifciiDA, a fluctuation with fc|| = also has uj = 0. If 
such modes are generated by the nonlinear interactions 
in the turbulence, this could be another possible cause for 
an excess of energy in their frequency spectra at u> ~ 0. 

Third, DM report that the frequency spectrum of their 
driving has the form P(uj) ~ 1/(oj 2 + u; 2 ), with ui c = 3 in 
their dimensionless units. The antenna power therefore 
scales P cx uj for uj <C uj c and P cx uj~ 2 for ui 3> ui c . How- 
ever, DM do not present a plot of the frequency spectrum 
of their driving, making it difficult to assess fairly the 
impact of the driving on the frequency spectra of their 



turbulence simulations. Since the parallel energy spec- 
trum of strong MHD turbulence 3 ^ is predicted to scale as 
i?(fe||) cx fc^ 2 , the linear relationship between the paral- 
lel wavenumbcr and frequency implies a frequency spec- 
trum E((j) cx uj~ 2 . In weak turbulence, the frequency 
spectrum would be steeper yet. Therefore, it is question- 
able whether one could observe even the strong turbulent 
frequency spectrum in the presence of their driving. 

One may attempt to judge the impact of the driving by 
examining Figure 8 in DM, a comparison of the frequency 
spectra from a driven and an undriven simulation. It is 
important to note that both are weak turbulence sim- 
ulations, so the contribution to the spectra due to the 
nonlinear turbulent fluctuations should be similar. The 
driven simulation shows a significantly larger signal over 
the frequency range 0.4 < uj < 10, suggesting that the 
driving has a significant, if not dominant, effect on the 
frequency spectrum over this range. DM attribute the 
broadband nature of the frequency spectra in their suite 
of simulations, presented in their Figure 2, to the inher- 
ently nonlinear nature of MHD turbulence. We suggest 
that a more careful evaluation of the impact of their driv- 
ing on the frequency spectrum is required to establish the 
validity of their conclusion. 

We believe the arguments outlined above raise serious 
questions about the conclusions that DM reach regard- 
ing the nature of MHD turbulence, in particular the claim 
that linear mode properties play little role in the turbu- 
lent evolution. The concept of critical balance implicitly 
assumes that linear wave modes do play a role in plasma 
turbulence. The numerical evidence presented here for 
critical balance in the kinetic Alfven wave cascade of dis- 
sipation range turbulence therefore indirectly supports 
the notion that linear wave modes do indeed play a role 
in strong plasma turbulence. 

VII. CONCLUSION 

We have developed a theoretical model for critically 
balanced Alfvcn/KAW turbulence and compared the re- 
sults of a fully nonlinear, driven, gyrokinetic simulation 
to the theoretical prediction. We find excellent qualita- 
tive and quantitative agreement with the predictions of 
critical balance in the dissipation range of KAW turbu- 
lence. This result constitutes the first evidence of critical 
balance to be observed in a damped and dissipative ki- 
netic turbulence simulation. 

Having found agreement with critical balance implies 
the anisotropic cascade of Alfven waves in the incr- 
tial range extends into the dissipation range, where the 
anisotropic scaling of KAW turbulence is observed to be 
approximately feii cx fc^j_ with < £ < 1/3. The upper 
bound of this result agrees with theoretical predictions 
for the dissipation range scaling based upon fluid descrip- 
tions, and the damping present in our kinetic simulation 
is expected to strengthen the anisotropy 3 ^. 

The results of our Eulerian frequency analysis per- 



9 



formed by temporally sampling the magnetic field data 
across a fixed x-y plane provide further evidence of a crit- 
ically balanced cascade of energy beginning at our driving 
frequency. Aside from the expected population of energy 
in the ui = mode, we find no evidence of excess energy 

pit.fipr hplow flip lowest or ahnvo the highest linpar mnrlps 
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Mon. Not. Roy. Astron. Soc, 407, L31-L35 (2010 
arXiv: 1002 .2096 [physics . spacc-ph] 
ld Q. Y. Luo and D. J. Wu, "Obscn 



3 Q. Y. Luo and D. J. Wu, "Observations of Anisotropic Scaling 
of Sola r Wind Turbulence," Astrophys. J. Lett., 714, L138-L141 



resolved in our simulation, in contradiction to Dmitruk 
and Matthaeus—. 

Although the current analysis does not directly ad- 
dress the importance linear wave theory to plasma turbu- 
lence, agreement with critical balance implies linear wave 
modes play some role in governing turbulence. Forth- 
coming simulation analyses will explore the importance 
of linear wave modes in fully developed, strong plasma 
turbulence in more detail. 
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